home *** CD-ROM | disk | FTP | other *** search
/ Aminet 30 / Aminet 30 (1999)(Schatztruhe)[!][Apr 1999].iso / Aminet / gfx / misc / gnuplot-3.7src.lha / gnuplot-3.7src / gnuplot-3.7.lha / gnuplot-3.7 / standard.c < prev    next >
C/C++ Source or Header  |  1998-12-08  |  25KB  |  1,103 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: standard.c,v 1.24 1998/04/14 00:16:21 drd Exp $";
  3. #endif
  4.  
  5. /* GNUPLOT - standard.c */
  6.  
  7. /*[
  8.  * Copyright 1986 - 1993, 1998   Thomas Williams, Colin Kelley
  9.  *
  10.  * Permission to use, copy, and distribute this software and its
  11.  * documentation for any purpose with or without fee is hereby granted,
  12.  * provided that the above copyright notice appear in all copies and
  13.  * that both that copyright notice and this permission notice appear
  14.  * in supporting documentation.
  15.  *
  16.  * Permission to modify the software is granted, but not the right to
  17.  * distribute the complete modified source code.  Modifications are to
  18.  * be distributed as patches to the released version.  Permission to
  19.  * distribute binaries produced by compiling modified sources is granted,
  20.  * provided you
  21.  *   1. distribute the corresponding source modifications from the
  22.  *    released version in the form of a patch file along with the binaries,
  23.  *   2. add special version identification to distinguish your version
  24.  *    in addition to the base release version number,
  25.  *   3. provide your name and address as the primary contact for the
  26.  *    support of your modified version, and
  27.  *   4. retain our contact information in regard to use of the base
  28.  *    software.
  29.  * Permission to distribute the released version of the source code along
  30.  * with corresponding source modifications in the form of a patch file is
  31.  * granted with same provisions 2 through 4 for binary distributions.
  32.  *
  33.  * This software is provided "as is" without express or implied warranty
  34.  * to the extent permitted by applicable law.
  35. ]*/
  36.  
  37. #include "plot.h"
  38. #include "setshow.h" /* for ang2rad */
  39. #include "fnproto.h"
  40.  
  41. extern struct value stack[STACK_DEPTH];
  42. extern int s_p;
  43. extern double zero;
  44.  
  45. static double jzero __PROTO((double x));
  46. static double pzero __PROTO((double x));
  47. static double qzero __PROTO((double x));
  48. static double yzero __PROTO((double x));
  49. static double rj0 __PROTO((double x));
  50. static double ry0 __PROTO((double x));
  51. static double jone __PROTO((double x));
  52. static double pone __PROTO((double x));
  53. static double qone __PROTO((double x));
  54. static double yone __PROTO((double x));
  55. static double rj1 __PROTO((double x));
  56. static double ry1 __PROTO((double x));
  57.  
  58. /* The bessel function approximations here are from
  59.  * "Computer Approximations"
  60.  * by Hart, Cheney et al.
  61.  * John Wiley & Sons, 1968
  62.  */
  63.  
  64. /* There appears to be a mistake in Hart, Cheney et al. on page 149.
  65.  * Where it list Qn(x)/x ~ P(z*z)/Q(z*z), z = 8/x, it should read
  66.  *               Qn(x)/z ~ P(z*z)/Q(z*z), z = 8/x
  67.  * In the functions below, Qn(x) is implementated using the later
  68.  * equation.
  69.  * These bessel functions are accurate to about 1e-13
  70.  */
  71.  
  72. #if (defined (ATARI) || defined (MTOS)) && defined(__PUREC__)
  73. /* Sorry. But PUREC bugs here.
  74.  * These bessel functions are NOT accurate to about 1e-13
  75.  */
  76.  
  77. #define PI_ON_FOUR     0.785398163397448309615661
  78. #define PI_ON_TWO     1.570796326794896619231313
  79. #define THREE_PI_ON_FOUR 2.356194490192344928846982
  80. #define TWO_ON_PI     0.636619772367581343075535
  81.  
  82. static double dzero = 0.0;
  83.  
  84. /* jzero for x in [0,8]
  85.  * Index 5849, 19.22 digits precision
  86.  */
  87. static double pjzero[] = {
  88.      0.493378725179413356181681e+21,
  89.     -0.117915762910761053603844e+21,
  90.      0.638205934107235656228943e+19,
  91.     -0.136762035308817138686542e+18,
  92.      0.143435493914034611166432e+16,
  93.     -0.808522203485379387119947e+13,
  94.      0.250715828553688194555516e+11,
  95.     -0.405041237183313270636066e+8,
  96.      0.268578685698001498141585e+5
  97. };
  98.  
  99. static double qjzero[] = {
  100.      0.493378725179413356211328e+21,
  101.      0.542891838409228516020019e+19,
  102.      0.302463561670946269862733e+17,
  103.      0.112775673967979850705603e+15,
  104.      0.312304311494121317257247e+12,
  105.      0.669998767298223967181403e+9,
  106.      0.111463609846298537818240e+7,
  107.      0.136306365232897060444281e+4,
  108.      0.1e+1
  109. };
  110.  
  111. /* pzero for x in [8,inf]
  112.  * Index 6548, 18.16 digits precision
  113.  */
  114. static double ppzero[] = {
  115.      0.227790901973046843022700e+5,
  116.      0.413453866395807657967802e+5,
  117.      0.211705233808649443219340e+5,
  118.      0.348064864432492703474453e+4,
  119.      0.153762019090083542957717e+3,
  120.      0.889615484242104552360748e+0
  121. };
  122.  
  123. static double qpzero[] = {
  124.      0.227790901973046843176842e+5,
  125.      0.413704124955104166398920e+5,
  126.      0.212153505618801157304226e+5,
  127.      0.350287351382356082073561e+4,
  128.      0.157111598580808936490685e+3,
  129.      0.1e+1
  130. };
  131.  
  132. /* qzero for x in [8,inf]
  133.  * Index 6948, 18.33 digits precision
  134.  */
  135. static double pqzero[] = {
  136.     -0.892266002008000940984692e+2,
  137.     -0.185919536443429938002522e+3,
  138.     -0.111834299204827376112621e+3,
  139.     -0.223002616662141984716992e+2,
  140.     -0.124410267458356384591379e+1,
  141.     -0.8803330304868075181663e-2,
  142. };
  143.  
  144. static double qqzero[] = {
  145.      0.571050241285120619052476e+4,
  146.      0.119511315434346136469526e+5,
  147.      0.726427801692110188369134e+4,
  148.      0.148872312322837565816135e+4,
  149.      0.905937695949931258588188e+2,
  150.      0.1e+1
  151. };
  152.  
  153.  
  154. /* yzero for x in [0,8]
  155.  * Index 6245, 18.78 digits precision
  156.  */
  157. static double pyzero[] = {
  158.     -0.275028667862910958370193e+20,
  159.      0.658747327571955492599940e+20,
  160.     -0.524706558111276494129735e+19,
  161.      0.137562431639934407857134e+18,
  162.     -0.164860581718572947312208e+16,
  163.      0.102552085968639428450917e+14,
  164.     -0.343637122297904037817103e+11,
  165.      0.591521346568688965427383e+8,
  166.     -0.413703549793314855412524e+5
  167. };
  168.  
  169. static double qyzero[] = {
  170.      0.372645883898616588198998e+21,
  171.      0.419241704341083997390477e+19,
  172.      0.239288304349978185743936e+17,
  173.      0.916203803407518526248915e+14,
  174.      0.261306575504108124956848e+12,
  175.      0.579512264070072953738009e+9,
  176.      0.100170264128890626566665e+7,
  177.      0.128245277247899380417633e+4,
  178.      0.1e+1
  179. };
  180.  
  181.  
  182. /* jone for x in [0,8]
  183.  * Index 6050, 20.98 digits precision
  184.  */
  185. static double pjone[] = {
  186.      0.581199354001606143928051e+21,
  187.     -0.667210656892491629802094e+20,
  188.      0.231643358063400229793182e+19,
  189.     -0.358881756991010605074364e+17,
  190.      0.290879526383477540973760e+15,
  191.     -0.132298348033212645312547e+13,
  192.      0.341323418230170053909129e+10,
  193.     -0.469575353064299585976716e+7,
  194.      0.270112271089232341485679e+4
  195. };
  196.  
  197. static double qjone[] = {
  198.      0.116239870800321228785853e+22,
  199.      0.118577071219032099983711e+20,
  200.      0.609206139891752174610520e+17,
  201.      0.208166122130760735124018e+15,
  202.      0.524371026216764971540673e+12,
  203.      0.101386351435867398996705e+10,
  204.      0.150179359499858550592110e+7,
  205.      0.160693157348148780197092e+4,
  206.      0.1e+1
  207. };
  208.  
  209.  
  210. /* pone for x in [8,inf]
  211.  * Index 6749, 18.11 digits precision
  212.  */
  213. static double ppone[] = {
  214.      0.352246649133679798341724e+5,
  215.      0.627588452471612812690057e+5,
  216.      0.313539631109159574238670e+5,
  217.      0.498548320605943384345005e+4,
  218.      0.211152918285396238210572e+3,
  219.      0.12571716929145341558495e+1
  220. };
  221.  
  222. static double qpone[] = {
  223.      0.352246649133679798068390e+5,
  224.      0.626943469593560511888834e+5,
  225.      0.312404063819041039923016e+5,
  226.      0.493039649018108897938610e+4,
  227.      0.203077518913475932229357e+3,
  228.      0.1e+1
  229. };
  230.  
  231. /* qone for x in [8,inf]
  232.  * Index 7149, 18.28 digits precision
  233.  */
  234. static double pqone[] = {
  235.      0.351175191430355282253332e+3,
  236.      0.721039180490447503928086e+3,
  237.      0.425987301165444238988699e+3,
  238.      0.831898957673850827325226e+2,
  239.      0.45681716295512267064405e+1,
  240.      0.3532840052740123642735e-1
  241. };
  242.  
  243. static double qqone[] = {
  244.      0.749173741718091277145195e+4,
  245.      0.154141773392650970499848e+5,
  246.      0.915223170151699227059047e+4,
  247.      0.181118670055235135067242e+4,
  248.      0.103818758546213372877664e+3,
  249.      0.1e+1
  250. };
  251.  
  252.  
  253. /* yone for x in [0,8]
  254.  * Index 6444, 18.24 digits precision
  255.  */
  256. static double pyone[] = {
  257.     -0.292382196153296254310105e+20,
  258.      0.774852068218683964508809e+19,
  259.     -0.344104806308411444618546e+18,
  260.      0.591516076049007061849632e+16,
  261.     -0.486331694256717507482813e+14,
  262.      0.204969667374566218261980e+12,
  263.     -0.428947196885524880182182e+9,
  264.      0.355692400983052605669132e+6
  265. };
  266.  
  267. static double qyone[] = {
  268.      0.149131151130292035017408e+21,
  269.      0.181866284170613498688507e+19,
  270.      0.113163938269888452690508e+17,
  271.      0.475517358888813771309277e+14,
  272.      0.150022169915670898716637e+12,
  273.      0.371666079862193028559693e+9,
  274.      0.726914730719888456980191e+6,
  275.      0.107269614377892552332213e+4,
  276.      0.1e+1
  277. };
  278.  
  279. #else
  280.  
  281. #define PI_ON_FOUR       0.78539816339744830961566084581987572
  282. #define PI_ON_TWO        1.57079632679489661923131269163975144
  283. #define THREE_PI_ON_FOUR 2.35619449019234492884698253745962716
  284. #define TWO_ON_PI        0.63661977236758134307553505349005744
  285.  
  286. static double dzero = 0.0;
  287.  
  288. /* jzero for x in [0,8]
  289.  * Index 5849, 19.22 digits precision
  290.  */
  291. static double GPFAR pjzero[] = {
  292.      0.4933787251794133561816813446e+21,
  293.     -0.11791576291076105360384408e+21,
  294.      0.6382059341072356562289432465e+19,
  295.     -0.1367620353088171386865416609e+18,
  296.      0.1434354939140346111664316553e+16,
  297.     -0.8085222034853793871199468171e+13,
  298.      0.2507158285536881945555156435e+11,
  299.     -0.4050412371833132706360663322e+8,
  300.      0.2685786856980014981415848441e+5
  301. };
  302.  
  303. static double GPFAR qjzero[] = {
  304.     0.4933787251794133562113278438e+21,
  305.     0.5428918384092285160200195092e+19,
  306.     0.3024635616709462698627330784e+17,
  307.     0.1127756739679798507056031594e+15,
  308.     0.3123043114941213172572469442e+12,
  309.     0.669998767298223967181402866e+9,
  310.     0.1114636098462985378182402543e+7,
  311.     0.1363063652328970604442810507e+4,
  312.     0.1e+1
  313. };
  314.  
  315. /* pzero for x in [8,inf]
  316.  * Index 6548, 18.16 digits precision
  317.  */
  318. static double GPFAR ppzero[] = {
  319.     0.2277909019730468430227002627e+5,
  320.     0.4134538663958076579678016384e+5,
  321.     0.2117052338086494432193395727e+5,
  322.     0.348064864432492703474453111e+4,
  323.     0.15376201909008354295771715e+3,
  324.     0.889615484242104552360748e+0
  325. };
  326.  
  327. static double GPFAR qpzero[] = {
  328.     0.2277909019730468431768423768e+5,
  329.     0.4137041249551041663989198384e+5,
  330.     0.2121535056188011573042256764e+5,
  331.     0.350287351382356082073561423e+4,
  332.     0.15711159858080893649068482e+3,
  333.     0.1e+1
  334. };
  335.  
  336. /* qzero for x in [8,inf]
  337.  * Index 6948, 18.33 digits precision
  338.  */
  339. static double GPFAR pqzero[] = {
  340.     -0.8922660020080009409846916e+2,
  341.     -0.18591953644342993800252169e+3,
  342.     -0.11183429920482737611262123e+3,
  343.     -0.2230026166621419847169915e+2,
  344.     -0.124410267458356384591379e+1,
  345.     -0.8803330304868075181663e-2,
  346. };
  347.  
  348. static double GPFAR qqzero[] = {
  349.     0.571050241285120619052476459e+4,
  350.     0.1195113154343461364695265329e+5,
  351.     0.726427801692110188369134506e+4,
  352.     0.148872312322837565816134698e+4,
  353.     0.9059376959499312585881878e+2,
  354.     0.1e+1
  355. };
  356.  
  357.  
  358. /* yzero for x in [0,8]
  359.  * Index 6245, 18.78 digits precision
  360.  */
  361. static double GPFAR pyzero[] = {
  362.     -0.2750286678629109583701933175e+20,
  363.      0.6587473275719554925999402049e+20,
  364.     -0.5247065581112764941297350814e+19,
  365.      0.1375624316399344078571335453e+18,
  366.     -0.1648605817185729473122082537e+16,
  367.      0.1025520859686394284509167421e+14,
  368.     -0.3436371222979040378171030138e+11,
  369.      0.5915213465686889654273830069e+8,
  370.     -0.4137035497933148554125235152e+5
  371. };
  372.  
  373. static double GPFAR qyzero[] = {
  374.     0.3726458838986165881989980739e+21,
  375.     0.4192417043410839973904769661e+19,
  376.     0.2392883043499781857439356652e+17,
  377.     0.9162038034075185262489147968e+14,
  378.     0.2613065755041081249568482092e+12,
  379.     0.5795122640700729537380087915e+9,
  380.     0.1001702641288906265666651753e+7,
  381.     0.1282452772478993804176329391e+4,
  382.     0.1e+1
  383. };
  384.  
  385.  
  386. /* jone for x in [0,8]
  387.  * Index 6050, 20.98 digits precision
  388.  */
  389. static double GPFAR pjone[] = {
  390.      0.581199354001606143928050809e+21,
  391.     -0.6672106568924916298020941484e+20,
  392.      0.2316433580634002297931815435e+19,
  393.     -0.3588817569910106050743641413e+17,
  394.      0.2908795263834775409737601689e+15,
  395.     -0.1322983480332126453125473247e+13,
  396.      0.3413234182301700539091292655e+10,
  397.     -0.4695753530642995859767162166e+7,
  398.      0.270112271089232341485679099e+4
  399. };
  400.  
  401. static double GPFAR qjone[] = {
  402.     0.11623987080032122878585294e+22,
  403.     0.1185770712190320999837113348e+20,
  404.     0.6092061398917521746105196863e+17,
  405.     0.2081661221307607351240184229e+15,
  406.     0.5243710262167649715406728642e+12,
  407.     0.1013863514358673989967045588e+10,
  408.     0.1501793594998585505921097578e+7,
  409.     0.1606931573481487801970916749e+4,
  410.     0.1e+1
  411. };
  412.  
  413.  
  414. /* pone for x in [8,inf]
  415.  * Index 6749, 18.11 digits precision
  416.  */
  417. static double GPFAR ppone[] = {
  418.     0.352246649133679798341724373e+5,
  419.     0.62758845247161281269005675e+5,
  420.     0.313539631109159574238669888e+5,
  421.     0.49854832060594338434500455e+4,
  422.     0.2111529182853962382105718e+3,
  423.     0.12571716929145341558495e+1
  424. };
  425.  
  426. static double GPFAR qpone[] = {
  427.     0.352246649133679798068390431e+5,
  428.     0.626943469593560511888833731e+5,
  429.     0.312404063819041039923015703e+5,
  430.     0.4930396490181088979386097e+4,
  431.     0.2030775189134759322293574e+3,
  432.     0.1e+1
  433. };
  434.  
  435. /* qone for x in [8,inf]
  436.  * Index 7149, 18.28 digits precision
  437.  */
  438. static double GPFAR pqone[] = {
  439.     0.3511751914303552822533318e+3,
  440.     0.7210391804904475039280863e+3,
  441.     0.4259873011654442389886993e+3,
  442.     0.831898957673850827325226e+2,
  443.     0.45681716295512267064405e+1,
  444.     0.3532840052740123642735e-1
  445. };
  446.  
  447. static double GPFAR qqone[] = {
  448.     0.74917374171809127714519505e+4,
  449.     0.154141773392650970499848051e+5,
  450.     0.91522317015169922705904727e+4,
  451.     0.18111867005523513506724158e+4,
  452.     0.1038187585462133728776636e+3,
  453.     0.1e+1
  454. };
  455.  
  456.  
  457. /* yone for x in [0,8]
  458.  * Index 6444, 18.24 digits precision
  459.  */
  460. static double GPFAR pyone[] = {
  461.     -0.2923821961532962543101048748e+20,
  462.      0.7748520682186839645088094202e+19,
  463.     -0.3441048063084114446185461344e+18,
  464.      0.5915160760490070618496315281e+16,
  465.     -0.4863316942567175074828129117e+14,
  466.      0.2049696673745662182619800495e+12,
  467.     -0.4289471968855248801821819588e+9,
  468.      0.3556924009830526056691325215e+6
  469. };
  470.  
  471. static double GPFAR qyone[] = {
  472.     0.1491311511302920350174081355e+21,
  473.     0.1818662841706134986885065935e+19,
  474.     0.113163938269888452690508283e+17,
  475.     0.4755173588888137713092774006e+14,
  476.     0.1500221699156708987166369115e+12,
  477.     0.3716660798621930285596927703e+9,
  478.     0.726914730719888456980191315e+6,
  479.     0.10726961437789255233221267e+4,
  480.     0.1e+1
  481. };
  482.  
  483. #endif /* (ATARI || MTOS) && __PUREC__ */
  484.  
  485. void f_real()
  486. {
  487. struct value a;
  488.     push( Gcomplex(&a,real(pop(&a)), 0.0) );
  489. }
  490.  
  491. void f_imag()
  492. {
  493. struct value a;
  494.     push( Gcomplex(&a,imag(pop(&a)), 0.0) );
  495. }
  496.  
  497.  
  498. /* ang2rad is 1 if we are in radians, or pi/180 if we are in degrees */
  499.  
  500. void f_arg()
  501. {
  502. struct value a;
  503.     push( Gcomplex(&a,angle(pop(&a))/ang2rad, 0.0) );
  504. }
  505.  
  506. void f_conjg()
  507. {
  508. struct value a;
  509.     (void) pop(&a);
  510.     push( Gcomplex(&a,real(&a),-imag(&a) ));
  511. }
  512.  
  513. /* ang2rad is 1 if we are in radians, or pi/180 if we are in degrees */
  514.  
  515. void f_sin()
  516. {
  517. struct value a;
  518.     (void) pop(&a);
  519.     push( Gcomplex(&a,sin(ang2rad*real(&a))*cosh(ang2rad*imag(&a)), cos(ang2rad*real(&a))*sinh(ang2rad*imag(&a))) );
  520. }
  521.  
  522. void f_cos()
  523. {
  524. struct value a;
  525.     (void) pop(&a);
  526.     push( Gcomplex(&a,cos(ang2rad*real(&a))*cosh(ang2rad*imag(&a)), -sin(ang2rad*real(&a))*sinh(ang2rad*imag(&a))));
  527. }
  528.  
  529. void f_tan()
  530. {
  531. struct value a;
  532. register double den;
  533.     (void) pop(&a);
  534.     if (imag(&a) == 0.0)
  535.         push( Gcomplex(&a,tan(ang2rad*real(&a)),0.0) );
  536.     else {
  537.         den = cos(2*ang2rad*real(&a))+cosh(2*ang2rad*imag(&a));
  538.         if (den == 0.0) {
  539.             undefined = TRUE;
  540.             push( &a );
  541.         }
  542.         else
  543.             push( Gcomplex(&a,sin(2*ang2rad*real(&a))/den, sinh(2*ang2rad*imag(&a))/den) );
  544.     }
  545. }
  546.  
  547. void f_asin()
  548. {
  549. struct value a;
  550. register double alpha, beta, x, y;
  551.     (void) pop(&a);
  552.     x = real(&a); y = imag(&a);
  553.     if (y == 0.0 && fabs(x) <= 1.0) {
  554.         push( Gcomplex(&a,asin(x)/ang2rad,0.0) );
  555.     } else if (x == 0.0) {
  556.         push( Gcomplex(&a, 0.0, -log(-y+sqrt(y*y+1))/ang2rad) );
  557.     } else {
  558.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  559.         if (beta > 1) beta = 1;    /* Avoid rounding error problems */
  560.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  561.         push( Gcomplex(&a,asin(beta)/ang2rad, -log(alpha + sqrt(alpha*alpha-1))/ang2rad) );
  562.     }
  563. }
  564.  
  565. void f_acos()
  566. {
  567. struct value a;
  568. register double x, y;
  569.     (void) pop(&a);
  570.     x = real(&a); y = imag(&a);
  571.     if (y == 0.0 && fabs(x) <= 1.0) {
  572.             /* real result */
  573.             push( Gcomplex(&a,acos(x)/ang2rad,0.0) );
  574.     } else { 
  575.         double alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  576.         double beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  577.         if (beta > 1)
  578.             beta = 1;    /* Avoid rounding error problems */
  579.         else if (beta < -1)
  580.             beta = -1;
  581.         push( Gcomplex(&a,acos(beta)/ang2rad, log(alpha + sqrt(alpha*alpha-1))/ang2rad));
  582.     }
  583. }
  584.  
  585. void f_atan()
  586. {
  587. struct value a;
  588. register double x, y, u, v, w, z;
  589.     (void) pop(&a);
  590.     x = real(&a); y = imag(&a);
  591.     if (y == 0.0)
  592.         push( Gcomplex(&a,atan(x)/ang2rad, 0.0) );
  593.     else if (x == 0.0 && fabs(y) >= 1.0) {
  594.         undefined = TRUE;
  595.         push(Gcomplex(&a,0.0, 0.0));
  596.     } else {
  597.             if (x >= 0) {
  598.                 u = x;
  599.             v = y;
  600.         } else {
  601.                 u = -x;
  602.             v = -y;
  603.         }
  604.         
  605.             z = atan(2*u/(1-u*u-v*v));
  606.         w = log((u*u+(v+1)*(v+1))/(u*u+(v-1)*(v-1)))/4;
  607.         if (z < 0)
  608.                 z = z + 2*PI_ON_TWO;
  609.         if (x < 0) {
  610.                 z = -z;
  611.             w = -w;
  612.         }
  613.         push( Gcomplex(&a,0.5*z/ang2rad, w) );
  614.     }
  615. }
  616.  
  617. /* real parts only */
  618. void f_atan2()
  619. {
  620.   struct value a;
  621.   double x;
  622.   double y;
  623.  
  624.   x = real(pop(&a));
  625.   y = real(pop(&a));
  626.  
  627.   if (x == 0.0 && y == 0.0) {
  628.     undefined = TRUE;
  629.     push(Ginteger(&a,0));
  630.   }
  631.   push(Gcomplex(&a,atan2(y,x)/ang2rad,0.0));
  632. }
  633.  
  634.  
  635. void f_sinh()
  636. {
  637. struct value a;
  638.     (void) pop(&a);
  639.     push( Gcomplex(&a,sinh(real(&a))*cos(imag(&a)), cosh(real(&a))*sin(imag(&a))) );
  640. }
  641.  
  642. void f_cosh()
  643. {
  644. struct value a;
  645.     (void) pop(&a);
  646.     push( Gcomplex(&a,cosh(real(&a))*cos(imag(&a)), sinh(real(&a))*sin(imag(&a))) );
  647. }
  648.  
  649. void f_tanh()
  650. {
  651. struct value a;
  652. register double den;
  653.     (void) pop(&a);
  654.     den = cosh(2*real(&a)) + cos(2*imag(&a));
  655.     push( Gcomplex(&a,sinh(2*real(&a))/den, sin(2*imag(&a))/den) );
  656. }
  657.  
  658. void f_asinh()
  659. {
  660. struct value a;                          /* asinh(z) = -I*asin(I*z) */
  661. register double alpha, beta, x, y;
  662.     (void) pop(&a);
  663.     x = -imag(&a); y = real(&a);
  664.     if (y == 0.0 && fabs(x) <= 1.0) {
  665.             push( Gcomplex(&a, 0.0, -asin(x)/ang2rad) );
  666.     } else if (y == 0.0) {
  667.             push( Gcomplex(&a, 0.0, 0.0) );
  668.         undefined = TRUE;
  669.     } else if (x == 0.0) {
  670.             push( Gcomplex(&a, log(y+sqrt(y*y+1))/ang2rad, 0.0) );
  671.     } else {
  672.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  673.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  674.         push( Gcomplex(&a, log(alpha + sqrt(alpha*alpha-1))/ang2rad, -asin(beta)/ang2rad) );
  675.     }
  676. }
  677.  
  678. void f_acosh()
  679. {
  680. struct value a;
  681. register double alpha, beta, x, y;        /* acosh(z) = I*acos(z) */
  682.     (void) pop(&a);
  683.     x = real(&a); y = imag(&a);
  684.     if (y == 0.0 && fabs(x) <= 1.0) {
  685.         push( Gcomplex(&a, 0.0, acos(x)/ang2rad) );
  686.     } else if (y == 0) {
  687.             push( Gcomplex(&a, log(x+sqrt(x*x-1))/ang2rad, 0.0) );
  688.     } else {
  689.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  690.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  691.         push( Gcomplex(&a, log(alpha + sqrt(alpha*alpha-1))/ang2rad, acos(beta)/ang2rad));
  692.     }
  693. }
  694.  
  695. void f_atanh()
  696. {
  697. struct value a;
  698. register double x, y, u, v, w, z;        /* atanh(z) = -I*atan(I*z) */
  699.     (void) pop(&a);
  700.     x = -imag(&a); y = real(&a);
  701.     if (y == 0.0)
  702.         push( Gcomplex(&a, 0.0, -atan(x)/ang2rad) );
  703.     else if (x == 0.0 && fabs(y) >= 1.0) {
  704.         undefined = TRUE;
  705.         push(Gcomplex(&a,0.0, 0.0));
  706.     } else {
  707.             if (x >= 0) {
  708.                 u = x;
  709.             v = y;
  710.         } else {
  711.                 u = -x;
  712.             v = -y;
  713.         }
  714.         
  715.             z = atan(2*u/(1-u*u-v*v));
  716.         w = log((u*u+(v+1)*(v+1))/(u*u+(v-1)*(v-1)))/4;
  717.         if (z < 0)
  718.                 z = z + 2*PI_ON_TWO;
  719.         if (x < 0) {
  720.                 z = -z;
  721.             w = -w;
  722.         }
  723.         push( Gcomplex(&a, w, -0.5*z/ang2rad) );
  724.     }
  725. }
  726.  
  727. void f_int()
  728. {
  729. struct value a;
  730.     push( Ginteger(&a,(int)real(pop(&a))) );
  731. }
  732.  
  733.  
  734. void f_abs()
  735. {
  736. struct value a;
  737.     (void) pop(&a);
  738.     switch (a.type) {
  739.         case INTGR:
  740.             push( Ginteger(&a,abs(a.v.int_val)) );            
  741.             break;
  742.         case CMPLX:
  743.             push( Gcomplex(&a,magnitude(&a), 0.0) );
  744.     }
  745. }
  746.  
  747. void f_sgn()
  748. {
  749. struct value a;
  750.     (void) pop(&a);
  751.     switch(a.type) {
  752.         case INTGR:
  753.             push( Ginteger(&a,(a.v.int_val > 0) ? 1 : 
  754.                     (a.v.int_val < 0) ? -1 : 0) );
  755.             break;
  756.         case CMPLX:
  757.             push( Ginteger(&a,(a.v.cmplx_val.real > 0.0) ? 1 : 
  758.                     (a.v.cmplx_val.real < 0.0) ? -1 : 0) );
  759.             break;
  760.     }
  761. }
  762.  
  763.  
  764. void f_sqrt()
  765. {
  766. struct value a;
  767. register double mag;
  768.     (void) pop(&a);
  769.     mag = sqrt(magnitude(&a));
  770.     if (imag(&a) == 0.0) {
  771.         if (real(&a) < 0.0)
  772.             push( Gcomplex(&a,0.0,mag) );
  773.         else
  774.             push( Gcomplex(&a,mag, 0.0) );
  775.     } else {
  776.         /* -pi < ang < pi, so real(sqrt(z)) >= 0 */
  777.         double ang = angle(&a) / 2.0;
  778.         push( Gcomplex(&a,mag*cos(ang), mag*sin(ang)) );
  779.     }
  780. }
  781.  
  782.  
  783. void f_exp()
  784. {
  785. struct value a;
  786. register double mag, ang;
  787.     (void) pop(&a);
  788.     mag = gp_exp(real(&a));
  789.     ang = imag(&a);
  790.     push( Gcomplex(&a,mag*cos(ang), mag*sin(ang)) );
  791. }
  792.  
  793.  
  794. void f_log10()
  795. {
  796. struct value a;
  797.     (void) pop(&a);
  798.     push( Gcomplex(&a,log(magnitude(&a))/M_LN10, angle(&a)/M_LN10) );
  799. }
  800.  
  801.  
  802. void f_log()
  803. {
  804. struct value a;
  805.     (void) pop(&a);
  806.     push( Gcomplex(&a,log(magnitude(&a)), angle(&a)) );
  807. }
  808.  
  809.  
  810. void f_floor()
  811. {
  812. struct value a;
  813.  
  814.     (void) pop(&a);
  815.     switch (a.type) {
  816.         case INTGR:
  817.             push( Ginteger(&a,(int)floor((double)a.v.int_val)));            
  818.             break;
  819.         case CMPLX:
  820.             push( Ginteger(&a,(int)floor(a.v.cmplx_val.real)));
  821.     }
  822. }
  823.  
  824.  
  825. void f_ceil()
  826. {
  827. struct value a;
  828.  
  829.     (void) pop(&a);
  830.     switch (a.type) {
  831.         case INTGR:
  832.             push( Ginteger(&a,(int)ceil((double)a.v.int_val)));            
  833.             break;
  834.         case CMPLX:
  835.             push( Ginteger(&a,(int)ceil(a.v.cmplx_val.real)));
  836.     }
  837. }
  838.  
  839. /* bessel function approximations */
  840. static double jzero(x)
  841. double x;
  842. {
  843. double p, q, x2;
  844. int n;
  845.  
  846.     x2 = x * x;
  847.     p = pjzero[8];
  848.     q = qjzero[8];
  849.     for (n=7; n>=0; n--) {
  850.         p = p*x2 + pjzero[n];
  851.         q = q*x2 + qjzero[n];
  852.     }
  853.     return(p/q);
  854. }
  855.  
  856. static double pzero(x)
  857. double x;
  858. {
  859. double p, q, z, z2;
  860. int n;
  861.  
  862.     z = 8.0 / x;
  863.     z2 = z * z;
  864.     p = ppzero[5];
  865.     q = qpzero[5];
  866.     for (n=4; n>=0; n--) {
  867.         p = p*z2 + ppzero[n];
  868.         q = q*z2 + qpzero[n];
  869.     }
  870.     return(p/q);
  871. }
  872.  
  873. static double qzero(x)
  874. double x;
  875. {
  876. double p, q, z, z2;
  877. int n;
  878.  
  879.     z = 8.0 / x;
  880.     z2 = z * z;
  881.     p = pqzero[5];
  882.     q = qqzero[5];
  883.     for (n=4; n>=0; n--) {
  884.         p = p*z2 + pqzero[n];
  885.         q = q*z2 + qqzero[n];
  886.     }
  887.     return(p/q);
  888. }
  889.  
  890. static double yzero(x)
  891. double x;
  892. {
  893. double p, q, x2;
  894. int n;
  895.  
  896.     x2 = x * x;
  897.     p = pyzero[8];
  898.     q = qyzero[8];
  899.     for (n=7; n>=0; n--) {
  900.         p = p*x2 + pyzero[n];
  901.         q = q*x2 + qyzero[n];
  902.     }
  903.     return(p/q);
  904. }
  905.  
  906. static double rj0(x)
  907. double x;
  908. {
  909.     if ( x <= 0.0 )
  910.         x = -x;
  911.     if ( x < 8.0 )
  912.         return(jzero(x));
  913.     else
  914.         return( sqrt(TWO_ON_PI/x) *
  915.             (pzero(x)*cos(x-PI_ON_FOUR) - 8.0/x*qzero(x)*sin(x-PI_ON_FOUR)) );
  916.  
  917. }
  918.  
  919. static double ry0(x)
  920. double x;
  921. {
  922.     if ( x < 0.0 )
  923.         return(dzero/dzero); /* error */
  924.     if ( x < 8.0 )
  925.         return( yzero(x) + TWO_ON_PI*rj0(x)*log(x) );
  926.     else
  927.         return( sqrt(TWO_ON_PI/x) *
  928.             (pzero(x)*sin(x-PI_ON_FOUR) + 
  929.             (8.0/x)*qzero(x)*cos(x-PI_ON_FOUR)) );
  930.  
  931. }
  932.  
  933.  
  934. static double jone(x)
  935. double x;
  936. {
  937. double p, q, x2;
  938. int n;
  939.  
  940.     x2 = x * x;
  941.     p = pjone[8];
  942.     q = qjone[8];
  943.     for (n=7; n>=0; n--) {
  944.         p = p*x2 + pjone[n];
  945.         q = q*x2 + qjone[n];
  946.     }
  947.     return(p/q);
  948. }
  949.  
  950. static double pone(x)
  951. double x;
  952. {
  953. double p, q, z, z2;
  954. int n;
  955.  
  956.     z = 8.0 / x;
  957.     z2 = z * z;
  958.     p = ppone[5];
  959.     q = qpone[5];
  960.     for (n=4; n>=0; n--) {
  961.         p = p*z2 + ppone[n];
  962.         q = q*z2 + qpone[n];
  963.     }
  964.     return(p/q);
  965. }
  966.  
  967. static double qone(x)
  968. double x;
  969. {
  970. double p, q, z, z2;
  971. int n;
  972.  
  973.     z = 8.0 / x;
  974.     z2 = z * z;
  975.     p = pqone[5];
  976.     q = qqone[5];
  977.     for (n=4; n>=0; n--) {
  978.         p = p*z2 + pqone[n];
  979.         q = q*z2 + qqone[n];
  980.     }
  981.     return(p/q);
  982. }
  983.  
  984. static double yone(x)
  985. double x;
  986. {
  987. double p, q, x2;
  988. int n;
  989.  
  990.     x2 = x * x;
  991.     p = 0.0;
  992.     q = qyone[8];
  993.     for (n=7; n>=0; n--) {
  994.         p = p*x2 + pyone[n];
  995.         q = q*x2 + qyone[n];
  996.     }
  997.     return(p/q);
  998. }
  999.  
  1000. static double rj1(x)
  1001. double x;
  1002. {
  1003. double v,w;
  1004.     v = x;
  1005.     if ( x < 0.0 )
  1006.         x = -x;
  1007.     if ( x < 8.0 )
  1008.         return(v*jone(x));
  1009.     else {
  1010.         w = sqrt(TWO_ON_PI/x) *
  1011.             (pone(x)*cos(x-THREE_PI_ON_FOUR) - 
  1012.                8.0/x*qone(x)*sin(x-THREE_PI_ON_FOUR)) ;
  1013.         if (v < 0.0)
  1014.             w = -w;
  1015.         return( w );
  1016.     }
  1017. }
  1018.  
  1019. static double ry1(x)
  1020. double x;
  1021. {
  1022.     if ( x <= 0.0 )
  1023.         return(dzero/dzero); /* error */
  1024.     if ( x < 8.0 )
  1025.         return( x*yone(x) + TWO_ON_PI*(rj1(x)*log(x) - 1.0/x) );
  1026.     else
  1027.         return( sqrt(TWO_ON_PI/x) *
  1028.             (pone(x)*sin(x-THREE_PI_ON_FOUR) + 
  1029.             (8.0/x)*qone(x)*cos(x-THREE_PI_ON_FOUR)) );
  1030. }
  1031.  
  1032.  
  1033. void f_besj0()    
  1034. {
  1035. struct value a;
  1036.     (void) pop(&a);
  1037.     if (fabs(imag(&a)) > zero)
  1038.         int_error("can only do bessel functions of reals",NO_CARET);
  1039.     push( Gcomplex(&a,rj0(real(&a)),0.0) );
  1040. }
  1041.  
  1042.  
  1043. void f_besj1()    
  1044. {
  1045. struct value a;
  1046.     (void) pop(&a);
  1047.     if (fabs(imag(&a)) > zero)
  1048.         int_error("can only do bessel functions of reals",NO_CARET);
  1049.     push( Gcomplex(&a,rj1(real(&a)),0.0) );
  1050. }
  1051.  
  1052.  
  1053. void f_besy0()    
  1054. {
  1055. struct value a;
  1056.     (void) pop(&a);
  1057.     if (fabs(imag(&a)) > zero)
  1058.         int_error("can only do bessel functions of reals",NO_CARET);
  1059.     if (real(&a) > 0.0)
  1060.         push( Gcomplex(&a,ry0(real(&a)),0.0) );
  1061.     else {
  1062.         push( Gcomplex(&a,0.0,0.0) );
  1063.         undefined = TRUE ;
  1064.     }
  1065. }
  1066.  
  1067.  
  1068. void f_besy1()    
  1069. {
  1070. struct value a;
  1071.     (void) pop(&a);
  1072.     if (fabs(imag(&a)) > zero)
  1073.         int_error("can only do bessel functions of reals",NO_CARET);
  1074.     if (real(&a) > 0.0)
  1075.         push( Gcomplex(&a,ry1(real(&a)),0.0) );
  1076.     else {
  1077.         push( Gcomplex(&a,0.0,0.0) );
  1078.         undefined = TRUE ;
  1079.     }
  1080. }
  1081.  
  1082.  
  1083. /* functions for accessing fields from tm structure, for time series
  1084.  * they are all the same, so define a macro
  1085.  */
  1086.  
  1087. #define TIMEFUNC(name, field) \
  1088. void name() { \
  1089.   struct value a;  struct tm tm;         \
  1090.   (void) pop(&a);  ggmtime(&tm, real(&a)); \
  1091.   push(Gcomplex(&a, (double)tm.field, 0.0));        \
  1092. }
  1093.  
  1094. TIMEFUNC(f_tmsec, tm_sec)
  1095. TIMEFUNC(f_tmmin, tm_min)
  1096. TIMEFUNC(f_tmhour, tm_hour)
  1097. TIMEFUNC(f_tmmday, tm_mday)
  1098. TIMEFUNC(f_tmmon, tm_mon)
  1099. TIMEFUNC(f_tmyear, tm_year)
  1100. TIMEFUNC(f_tmwday, tm_wday)
  1101. TIMEFUNC(f_tmyday, tm_yday)
  1102.  
  1103.